www.gusucode.com > matlab编程NSCT分解 图像融合 各个融合指标评价体系 分解源码程序 > matlab编程NSCT分解 图像融合 各个融合指标评价体系 分解源码程序/NSCT/phasecongmono.m
% PHASECONGMONO - phase congruency of an image using monogenic filters
%
% This code is considerably faster than PHASECONG3 but you may prefer the
% output from PHASECONG3's oriented filters.
%
% There are potentially many arguments, here is the full usage:
%
% [PC or ft T] = ...
% phasecongmono(im, nscale, minWaveLength, mult, ...
% sigmaOnf, k, cutOff, g, deviationGain, noiseMethod)
%
% However, apart from the image, all parameters have defaults and the
% usage can be as simple as:
%
% phaseCong = phasecongmono(im);
%
% Arguments:
% Default values Description
%
% nscale 4 - Number of wavelet scales, try values 3-6
% A lower value will reveal more fine scale
% features. A larger value will highlight 'major'
% features.
% minWaveLength 3 - Wavelength of smallest scale filter.
% mult 2.1 - Scaling factor between successive filters.
% sigmaOnf 0.55 - Ratio of the standard deviation of the Gaussian
% describing the log Gabor filter's transfer function
% in the frequency domain to the filter center frequency.
% k 3.0 - No of standard deviations of the noise energy beyond
% the mean at which we set the noise threshold point.
% You may want to vary this up to a value of 10 or
% 20 for noisy images
% cutOff 0.5 - The fractional measure of frequency spread
% below which phase congruency values get penalized.
% g 10 - Controls the sharpness of the transition in
% the sigmoid function used to weight phase
% congruency for frequency spread.
% deviationGain 1.5 - Amplification to apply to the calculated phase
% deviation result. Increasing this sharpens the
% edge responses, but can also attenuate their
% magnitude if the gain is too large. Sensible
% values to use lie in the range 1-2.
% noiseMethod -1 - Parameter specifies method used to determine
% noise statistics.
% -1 use median of smallest scale filter responses
% -2 use mode of smallest scale filter responses
% 0+ use noiseMethod value as the fixed noise threshold
% A value of 0 will turn off all noise compensation.
%
% Returned values:
% PC - Phase congruency indicating edge significance
% or - Orientation image in integer degrees 0-180,
% positive anticlockwise.
% 0 corresponds to a vertical edge, 90 is horizontal.
% ft - Local weighted mean phase angle at every point in the
% image. A value of pi/2 corresponds to a bright line, 0
% corresponds to a step and -pi/2 is a dark line.
% T - Calculated noise threshold (can be useful for
% diagnosing noise characteristics of images). Once you know
% this you can then specify fixed thresholds and save some
% computation time.
%
% Notes on specifying parameters:
%
% The parameters can be specified as a full list eg.
% >> PC = phasecongmono(im, 5, 3, 2.5, 0.55, 2.0);
%
% or as a partial list with unspecified parameters taking on default values
% >> PC = phasecongmono(im, 5, 3);
%
% or as a partial list of parameters followed by some parameters specified via a
% keyword-value pair, remaining parameters are set to defaults, for example:
% >> PC = phasecongmono(im, 5, 3, 'k', 2.5);
%
% The convolutions are done via the FFT. Many of the parameters relate to the
% specification of the filters in the frequency plane. The values do not seem
% to be very critical and the defaults are usually fine. You may want to
% experiment with the values of 'nscales' and 'k', the noise compensation
% factor.
%
% Typical sequence of operations to obtain an edge image:
%
% >> [PC, or] = phasecongmono(imread('lena.tif'));
% >> nm = nonmaxsup(PC, or, 1.5); % nonmaxima suppression
% >> bw = hysthresh(nm, 0.1, 0.3); % hysteresis thresholding 0.1 - 0.3
% >> show(bw)
%
% Notes on filter settings to obtain even coverage of the spectrum
% sigmaOnf .85 mult 1.3
% sigmaOnf .75 mult 1.6 (filter bandwidth ~1 octave)
% sigmaOnf .65 mult 2.1
% sigmaOnf .55 mult 3 (filter bandwidth ~2 octaves)
%
% Note that better results are achieved using the large bandwidth filters.
% I generally use a sigmaOnf value of 0.55 or even smaller.
%
% See Also: PHASECONG, PHASECONG3, PHASESYMMONO, GABORCONVOLVE,
% PLOTGABORFILTERS, FILTERGRID
% References:
%
% Peter Kovesi, "Image Features From Phase Congruency". Videre: A
% Journal of Computer Vision Research. MIT Press. Volume 1, Number 3,
% Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html
%
% Michael Felsberg and Gerald Sommer, "A New Extension of Linear Signal
% Processing for Estimating Local Properties and Detecting Features". DAGM
% Symposium 2000, Kiel
%
% Michael Felsberg and Gerald Sommer. "The Monogenic Signal" IEEE
% Transactions on Signal Processing, 49(12):3136-3144, December 2001
%
% Peter Kovesi, "Phase Congruency Detects Corners and Edges". Proceedings
% DICTA 2003, Sydney Dec 10-12
% August 2008 Initial version developed from phasesymmono and phasecong2
% where local phase information is calculated via Monogenic
% filters. Simplification of noise compensation to speedup
% execution. Options to calculate noise statistics via median
% or mode of smallest filter response.
% April 2009 Return of T for 'instrumentation' of noise
% detection/compensation. Option to use a fixed threshold.
% Frequency width measure slightly improved.
% June 2009 20% Speed improvement through calculating phase deviation via
% acos() rather than computing cos(theta)-|sin(theta)| via dot
% and cross products. Phase congruency is formed properly in
% 2D rather than as multiple 1D computations as in
% phasecong3. Also, much smaller memory footprint.
% May 2013 Some tidying up, corrections to defualt parameters, changes
% to reflect my latest thinking of the final phase congruency
% calculation and addition of phase deviation gain parameter to
% sharpen up the output. Use of periodic fft.
% Copyright (c) 1996-2013 Peter Kovesi
% Centre for Exploration Targeting
% The University of Western Australia
% peter.kovesi at uwa edu au
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
function [PC, or, ft, T] = phasecongmono(varargin)
% Get arguments and/or default values
[im, nscale, minWaveLength, mult, sigmaOnf, k, ...
noiseMethod, cutOff, g, deviationGain] = checkargs(varargin(:));
epsilon = .0001; % Used to prevent division by zero.
[rows,cols] = size(im);
IM = perfft2(im); % Periodic Fourier transform of image
sumAn = zeros(rows,cols); % Matrix for accumulating filter response
% amplitude values.
sumf = zeros(rows,cols);
sumh1 = zeros(rows,cols);
sumh2 = zeros(rows,cols);
% Generate grid data for constructing filters in the frequency domain
[radius, u1, u2] = filtergrid(rows, cols);
% Get rid of the 0 radius value in the middle (at top left corner after
% fftshifting) so that taking the log of the radius, or dividing by the
% radius, will not cause trouble.
radius(1,1) = 1;
% Construct the monogenic filters in the frequency domain. The two
% filters would normally be constructed as follows
% H1 = i*u1./radius;
% H2 = i*u2./radius;
% However the two filters can be packed together as a complex valued
% matrix, one in the real part and one in the imaginary part. Do this by
% multiplying H2 by i and then adding it to H1 (note the subtraction
% because i*i = -1). When the convolution is performed via the fft the
% real part of the result will correspond to the convolution with H1 and
% the imaginary part with H2. This allows the two convolutions to be
% done as one in the frequency domain, saving time and memory.
H = (1i*u1 - u2)./radius;
% The two monogenic filters H1 and H2 are not selective in terms of the
% magnitudes of the frequencies. The code below generates bandpass
% log-Gabor filters which are point-wise multiplied by IM to produce
% different bandpass versions of the image before being convolved with H1
% and H2
% First construct a low-pass filter that is as large as possible, yet falls
% away to zero at the boundaries. All filters are multiplied by
% this to ensure no extra frequencies at the 'corners' of the FFT are
% incorporated as this can upset the normalisation process when
% calculating phase congruency
lp = lowpassfilter([rows,cols],.45,15); % Radius .4, 'sharpness' 15
for s = 1:nscale
wavelength = minWaveLength*mult^(s-1);
fo = 1.0/wavelength; % Centre frequency of filter.
logGabor = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));
logGabor = logGabor.*lp; % Apply low-pass filter
logGabor(1,1) = 0; % Set the value at the 0 frequency point of the
% filter back to zero (undo the radius fudge).
IMF = IM.*logGabor; % Bandpassed image in the frequency domain.
f = real(ifft2(IMF)); % Bandpassed image in spatial domain.
h = ifft2(IMF.*H); % Bandpassed monogenic filtering, real part of h contains
% convolution result with h1, imaginary part
% contains convolution result with h2.
h1 = real(h);
h2 = imag(h);
An = sqrt(f.^2 + h1.^2 + h2.^2); % Amplitude of this scale component.
sumAn = sumAn + An; % Sum of component amplitudes over scale.
sumf = sumf + f;
sumh1 = sumh1 + h1;
sumh2 = sumh2 + h2;
% At the smallest scale estimate noise characteristics from the
% distribution of the filter amplitude responses stored in sumAn.
% tau is the Rayleigh parameter that is used to describe the
% distribution.
if s == 1
if noiseMethod == -1 % Use median to estimate noise statistics
tau = median(sumAn(:))/sqrt(log(4));
elseif noiseMethod == -2 % Use mode to estimate noise statistics
tau = rayleighmode(sumAn(:));
end
maxAn = An;
else
% Record maximum amplitude of components across scales. This is needed
% to determine the frequency spread weighting.
maxAn = max(maxAn,An);
end
end % For each scale
% Form weighting that penalizes frequency distributions that are
% particularly narrow. Calculate fractional 'width' of the frequencies
% present by taking the sum of the filter response amplitudes and dividing
% by the maximum component amplitude at each point on the image. If
% there is only one non-zero component width takes on a value of 0, if
% all components are equal width is 1.
width = (sumAn./(maxAn + epsilon) - 1) / (nscale-1);
% Now calculate the sigmoidal weighting function.
weight = 1.0 ./ (1 + exp( (cutOff - width)*g));
% Automatically determine noise threshold
%
% Assuming the noise is Gaussian the response of the filters to noise will
% form Rayleigh distribution. We use the filter responses at the smallest
% scale as a guide to the underlying noise level because the smallest scale
% filters spend most of their time responding to noise, and only
% occasionally responding to features. Either the median, or the mode, of
% the distribution of filter responses can be used as a robust statistic to
% estimate the distribution mean and standard deviation as these are related
% to the median or mode by fixed constants. The response of the larger
% scale filters to noise can then be estimated from the smallest scale
% filter response according to their relative bandwidths.
%
% This code assumes that the expected reponse to noise on the phase
% congruency calculation is simply the sum of the expected noise responses
% of each of the filters. This is a simplistic overestimate, however these
% two quantities should be related by some constant that will depend on the
% filter bank being used. Appropriate tuning of the parameter 'k' will
% allow you to produce the desired output. (though the value of k seems to
% be not at all critical)
if noiseMethod >= 0 % We are using a fixed noise threshold
T = noiseMethod; % use supplied noiseMethod value as the threshold
else
% Estimate the effect of noise on the sum of the filter responses as
% the sum of estimated individual responses (this is a simplistic
% overestimate). As the estimated noise response at succesive scales
% is scaled inversely proportional to bandwidth we have a simple
% geometric sum.
totalTau = tau * (1 - (1/mult)^nscale)/(1-(1/mult));
% Calculate mean and std dev from tau using fixed relationship
% between these parameters and tau. See
% http://mathworld.wolfram.com/RayleighDistribution.html
EstNoiseEnergyMean = totalTau*sqrt(pi/2); % Expected mean and std
EstNoiseEnergySigma = totalTau*sqrt((4-pi)/2); % values of noise energy
T = EstNoiseEnergyMean + k*EstNoiseEnergySigma; % Noise threshold
end
%------ Final computation of key quantities -------
or = atan(-sumh2./sumh1); % Orientation - this varies +/- pi/2
or(or<0) = or(or<0)+pi; % Add pi to -ve orientation value so that all
% orientation values now range 0 - pi
or = fix(or/pi*180); % Quantize to 0 - 180 degrees (for NONMAXSUP)
ft = atan2(sumf,sqrt(sumh1.^2+sumh2.^2)); % Feature type - a phase angle
% -pi/2 to pi/2.
energy = sqrt(sumf.^2 + sumh1.^2 + sumh2.^2); % Overall energy
% Compute phase congruency. The original measure,
% PC = energy/sumAn
% is proportional to the weighted cos(phasedeviation). This is not very
% localised so this was modified to
% PC = cos(phasedeviation) - |sin(phasedeviation)|
% (Note this was actually calculated via dot and cross products.) This measure
% approximates
% PC = 1 - phasedeviation.
% However, rather than use dot and cross products it is simpler and more
% efficient to simply use acos(energy/sumAn) to obtain the weighted phase
% deviation directly. Note, in the expression below the noise threshold is
% not subtracted from energy immediately as this would interfere with the
% phase deviation computation. Instead it is applied as a weighting as a
% fraction by which energy exceeds the noise threshold. This weighting is
% applied in addition to the weighting for frequency spread. Note also the
% phase deviation gain factor which acts to sharpen up the edge response. A
% value of 1.5 seems to work well. Sensible values are from 1 to about 2.
PC = weight.*max(1 - deviationGain*acos(energy./(sumAn + epsilon)),0) ...
.* max(energy-T,0)./(energy+epsilon);
%------------------------------------------------------------------
% CHECKARGS
%
% Function to process the arguments that have been supplied, assign
% default values as needed and perform basic checks.
function [im, nscale, minWaveLength, mult, sigmaOnf, ...
k, noiseMethod, cutOff, g, deviationGain] = checkargs(arg)
nargs = length(arg);
if nargs < 1
error('No image supplied as an argument');
end
% Set up default values for all arguments and then overwrite them
% with with any new values that may be supplied
im = [];
nscale = 4; % Number of wavelet scales.
minWaveLength = 3; % Wavelength of smallest scale filter.
mult = 2.1; % Scaling factor between successive filters.
sigmaOnf = 0.55; % Ratio of the standard deviation of the
% Gaussian describing the log Gabor filter's
% transfer function in the frequency domain
% to the filter center frequency.
k = 3.0; % No of standard deviations of the noise
% energy beyond the mean at which we set the
% noise threshold point.
noiseMethod = -1; % Use the median response of smallest scale
% filter to estimate noise statistics
cutOff = 0.5;
g = 10;
deviationGain = 1.5;
% Allowed argument reading states
allnumeric = 1; % Numeric argument values in predefined order
keywordvalue = 2; % Arguments in the form of string keyword
% followed by numeric value
readstate = allnumeric; % Start in the allnumeric state
if readstate == allnumeric
for n = 1:nargs
if isa(arg{n},'char')
readstate = keywordvalue;
break;
else
if n == 1, im = arg{n};
elseif n == 2, nscale = arg{n};
elseif n == 3, minWaveLength = arg{n};
elseif n == 4, mult = arg{n};
elseif n == 5, sigmaOnf = arg{n};
elseif n == 6, k = arg{n};
elseif n == 7, cutOff = arg{n};
elseif n == 8, g = arg{n};
elseif n == 9, deviationGain = arg{n};
elseif n == 10, noiseMethod = arg{n};
end
end
end
end
% Code to handle parameter name - value pairs
if readstate == keywordvalue
while n < nargs
if ~isa(arg{n},'char') || ~isa(arg{n+1}, 'double')
error('There should be a parameter name - value pair');
end
if strncmpi(arg{n},'im' ,2), im = arg{n+1};
elseif strncmpi(arg{n},'nscale' ,2), nscale = arg{n+1};
elseif strncmpi(arg{n},'minWaveLength',2), minWaveLength = arg{n+1};
elseif strncmpi(arg{n},'mult' ,2), mult = arg{n+1};
elseif strncmpi(arg{n},'sigmaOnf',2), sigmaOnf = arg{n+1};
elseif strncmpi(arg{n},'k' ,1), k = arg{n+1};
elseif strncmpi(arg{n},'cutOff' ,2), cutOff = arg{n+1};
elseif strncmpi(arg{n},'g' ,1), g = arg{n+1};
elseif strncmpi(arg{n},'deviation',3), deviationGain = arg{n+1};
elseif strncmpi(arg{n},'noisemethod',3), noiseMethod = arg{n+1};
else error('Unrecognised parameter name');
end
n = n+2;
if n == nargs
error('Unmatched parameter name - value pair');
end
end
end
if isempty(im)
error('No image argument supplied');
end
if ndims(im) == 3
warning('Colour image supplied: converting image to greyscale...')
im = double(rgb2gray(im));
end
if ~isa(im, 'double')
im = double(im);
end
if nscale < 1
error('nscale must be an integer >= 1');
end
if minWaveLength < 2
error('It makes little sense to have a wavelength < 2');
end
%-------------------------------------------------------------------------
% RAYLEIGHMODE
%
% Computes mode of a vector/matrix of data that is assumed to come from a
% Rayleigh distribution.
%
% Usage: rmode = rayleighmode(data, nbins)
%
% Arguments: data - data assumed to come from a Rayleigh distribution
% nbins - Optional number of bins to use when forming histogram
% of the data to determine the mode.
%
% Mode is computed by forming a histogram of the data over 50 bins and then
% finding the maximum value in the histogram. Mean and standard deviation
% can then be calculated from the mode as they are related by fixed
% constants.
%
% mean = mode * sqrt(pi/2)
% std dev = mode * sqrt((4-pi)/2)
%
% See
% http://mathworld.wolfram.com/RayleighDistribution.html
% http://en.wikipedia.org/wiki/Rayleigh_distribution
%
function rmode = rayleighmode(data, nbins)
if nargin == 1
nbins = 50; % Default number of histogram bins to use
end
mx = max(data(:));
edges = 0:mx/nbins:mx;
n = histc(data(:),edges);
[dum,ind] = max(n); % Find maximum and index of maximum in histogram
rmode = (edges(ind)+edges(ind+1))/2;